home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1996 July / C-C++ Users Group Library July 1996.iso / vol_200 / 228_01 / mat.c < prev    next >
Encoding:
Text File  |  1987-07-29  |  6.3 KB  |  227 lines
  1. /*
  2. HEADER:         CUGXXX;
  3. TITLE:          Matrix decomposition & double-back substitution;
  4. DATE:           3-20-86;
  5. DESCRIPTION:    Matrix mathematics;
  6. KEYWORDS:       Matrix decomposition, Matrix double-back substitution;
  7. FILENAME:       MAT.C;
  8. WARNINGS:       None;
  9. AUTHORS:        Unknown;
  10. COMPILER:       DeSmet C;
  11. REFERENCES:     US-DISK 1308;
  12. ENDREF
  13. */
  14. mdecomp(a,n,intr1,cdim,det)  /* matrix decomposition routine
  15.                 in single precision */
  16.  
  17. /* includes partial pivoting and determinant calculation */
  18.  
  19. /* returns:    +1    if successful
  20.         -val    if singularity is detected, val is the step it was on
  21.         0    if input is invalid.
  22. */
  23. float a[]; /* a is the matrix stored in a[row][column] form */
  24. int n; /* n is the order of the matrix */
  25. int intr1[]; /* intr1 is a vector used to store the row interchanges */
  26. int cdim; /* cdim is the column dimension of a and intr1 */
  27. double *det; /* det is the determinant of the matrix */
  28.  
  29. {
  30.     int row,col,nm1; /* nm1 is n-1 */
  31.     int lpr; /* lpr is the largest pivot row */
  32.     int col2; /* col2 is the column for row interchanging */
  33.     float lpe; /* lpe is the largest pivot element */
  34.     float ab; /* ab is a tempary absolute value */
  35.     float temp; /* temp is a temporary storage */
  36.     double q; /* q is a quotient */
  37.     float *pa1,*pa2,*pa3,*pa4; /* 1st 2nd 3rd and 4th pointers to a */
  38.     int itr=1; /* itr keeps track of even or odd row interchanges */
  39.     int cdimp1; /* cdimp1 is cdim+1 */
  40.     double fabs();
  41.  
  42.     intr1[0] = 0; /* initialize intr1[0] to 0 to check errors in mback */
  43.     nm1 = n-1;
  44.     cdimp1=cdim+1;
  45.     *det = 0.;
  46.  
  47. /* check input */
  48.     if( n <= 1  ||    cdim < n )
  49.         return(0);
  50.  
  51. /* decompose matrix */
  52.     pa1 = a; /* pa1 points to a(col,col) */
  53.     for (col = 1 ; col <= nm1 ; col++)
  54.     {
  55.  
  56.         intr1[col] = 0;
  57.         lpr = col; /* the current row is col */
  58.         pa2 = pa1; /* pa2 points to a(row,col) */
  59.         lpe = fabs( *pa2 );
  60.         for (row = col+1 ; row <= n ; row++) /* search for largest pivot */
  61.         {
  62.         pa2 += cdim; /* pa2 points to a(row,col) */
  63.         ab = fabs( *pa2 );
  64.         /* if the largest pivot element is smaller than the current
  65.         row, then make the current row the pivot row */
  66.         if(lpe < ab)
  67.         {
  68.             lpr = row;
  69.             lpe = ab;
  70.         }
  71.         }
  72.  
  73.         if (lpe == 0.) /* check for singularity */
  74.         return(-col);
  75.  
  76.         if(lpr != col) /* check for necessity of row interchanges */
  77.         {
  78.         itr *= -1; /* each row interchange changes the sign of det */
  79.         intr1[col] = lpr;
  80.         pa2 = a + (lpr-1)*cdim; /* pa2 points to a(lpr,col2) */
  81.         pa3 = pa1 - col + 1; /* pa3 points to a(col,col2) */
  82.         for(col2 = 1 ; col2 <= n ; col2++) /* interchange rows */
  83.         {
  84.             temp = *pa2;
  85.             *pa2++ = *pa3;
  86.             *pa3++ = temp;
  87.         }
  88.         }
  89.  
  90.         pa3 = pa1; /* pa3 points to a(col,col2) */
  91.         for (col2 = col + 1 ; col2 <= n ; col2++) /* normalize */
  92.         *++pa3 /= *pa1;    /* a(col,col+1) through a(col,n) */
  93.  
  94.         pa2 = pa1; /* pa2 points to a(row,col) */
  95.         for (row = col+1 ; row <= n ; row++) /* calc new mat elements */
  96.         {
  97.         pa2 += cdim; /* pa2 points to a(row,col) */
  98.         pa3 = pa1; /* pa3 points to a(col,col2) */
  99.         pa4 = pa2; /* pa4 points to a(row,col2) */
  100.         for (col2 = col+1 ; col2 <= n ; col2++)
  101.             *++pa4 -= (*pa2)*(*++pa3);
  102.             /* the previous statement is a(row,col2) =
  103.             a(row,col2) - a(row,col)*a(col,col2) */
  104.         }
  105.         pa1 += cdimp1; /* pa1 points to a(col,col) */
  106.     }
  107.     if(*pa1 == 0.) /* check if a(n,n) is 0. */
  108.         return(-n);
  109.  
  110.  /* find determinant */
  111.  
  112.     *det = 1.;
  113.     pa1=a; /* pa1 points to a(col,col) */
  114.     for (col = 1 ; col <= n ; col++)
  115.     {
  116.         *det *= *pa1;
  117.         pa1 += cdimp1;
  118.     }
  119.     *det *= itr;
  120.     intr1[0] = itr;
  121.  
  122. /* return */
  123.     return(1);
  124.  
  125. }
  126.  
  127.  
  128.  
  129. mback(b,m,cdimb,intr1,a,n,cdima) /* matrix double back substitution
  130.                     in single precision */
  131. /* used with mdecomp */
  132. /* returns:
  133.         0    if input is invalid
  134.         +1    if successful
  135. */
  136.  
  137. float b[]; /* right hand side matrix on input, solution matrix on output */
  138. int m; /* number of right hand sides */
  139. int cdimb; /* column dimension of b */
  140. int intr1[]; /* vector used to store row interchanges */
  141. float a[]; /* decomposed a matrix, is not altered in this routine */
  142. int n; /* order of matrix a */
  143. int cdima; /* column dimension of a */
  144.  
  145. {
  146.     int row,col,bcol;
  147.     int nm1,bcolm1; /* nm1 is n-1, bcolm1 is bcol-1 */
  148.     int rowint; /* rowint stands for row interchange */
  149.     float *pa1,*pa2; /* pa1 and pa2 are pointers for a */
  150.     float *pb1,*pb2; /* pb1 and pb2 are pointers for b */
  151.     float temp; /* temporary storage */
  152.     int cdimap1; /* cdimap1 is cdimb+1 */
  153.  
  154. /* check input */
  155.     if (m < 1 || cdimb < m || n <= 1 || cdima < n || intr1[0] == 0)
  156.         return(0); /* intr1[0]=0 if there was an error in mdecomp or
  157.             a was singular */
  158.  
  159.     nm1 = n-1;
  160.  
  161. /* interchange the right hand side vector if row interchanges were made */
  162.  
  163.     for (row = 1 ; row <= nm1 ; row++)
  164.     {
  165.         rowint = intr1[row];
  166.         if (rowint != 0) /* there was a row interchange */
  167.         {
  168.         pb1 = b + (row-1)*cdimb; /* pb1 points to b(row,bcol) */
  169.         pb2 = b + (rowint-1)*cdimb; /* pb2 points to b(rowint,bcol) */
  170.         for (bcol = 1 ; bcol <= m ; bcol++)
  171.         {
  172.             temp = *pb1;
  173.             *pb1++ = *pb2;
  174.             *pb2++ = temp;
  175.         }
  176.         }
  177.     }
  178.  
  179. /* calculate new b's based on row addition
  180. as previously performed on a in mdecomp */
  181.  
  182.     cdimap1 = cdima + 1;
  183.  
  184.     for (bcol = 1 ; bcol <= m ; bcol++)
  185.     {
  186.         bcolm1 = bcol - 1;
  187.         pb1 = b + bcolm1; /* pb1 points to b(col,bcol) */
  188.         pa1 = a; /* pa1 points to a(col,col) */
  189.         for (col = 1 ; col <= n ; col++)
  190.         {
  191.         *pb1 /= *pa1;
  192.         pb2 = pb1; /* pb2 points to b(row,bcol) */
  193.         pa2 = pa1; /* pa2 points to a(row,col) */
  194.         for (row = col + 1 ; row <= n ; row++)
  195.         {
  196.             pb2 += cdimb;
  197.             pa2 += cdima;
  198.             *pb2 -= (*pa2)*(*pb1);
  199.         }
  200.         pb1 += cdimb;
  201.         pa1 += cdimap1;
  202.         }
  203.     }
  204.  
  205. /* perform the double back substitution */
  206.  
  207.     for (bcol = 1 ; bcol <= m ; bcol++)
  208.     {
  209.         bcolm1 = bcol - 1;
  210.         pb1 = b + n*cdimb + bcolm1; /* pb1 points to b(col+1,bcol) */
  211.         for (col = nm1 ; col >= 1 ; col--)
  212.         {
  213.         pa1 = a + col; /* pa1 points to a(row,col+1) */
  214.         pb1 -= cdimb;
  215.         pb2 = b + bcolm1; /* pb2 points to b(row,bcol) */
  216.         for (row = 1 ; row <= col ; row++)
  217.         {
  218.             *pb2 -= (*pa1)*(*pb1);
  219.             pa1 += cdima;
  220.             pb2 += cdimb;
  221.         }
  222.         }
  223.     }
  224.     return(1);
  225.  
  226. }
  227.